function [Type0,X_fixed,X_oo,X_rc,X_interact_i,X_interact_j,X_BVPdiscount_ij,X_dist_ij,typ_mat] = MCMCsetupMatrices(MCMC,Data)
%setup some matrices for MCMC
Data.nStudents = size(Data.Scores,1);
Data.nOptions = size(Data.MajorCode,1);

%% setup and preliminaries
Locj = Data.LocationXj;
Sinorm = Data.Sinorm;

%pricing stuff here:
XjArancel = Data.ArancelXj;   %listprice
AidAmount = Data.AidAmount;
Over600=sum(Data.Scores,2)/2>=3;
EverBVPj = Data.EverBVP;

Xi = Data.StudentXs(:, 1:2);
Xi(isnan(Xi)) = 0;
Loci = Data.LocationXi;
Xjnorm=Data.OptionXjnorm;
XjSelectivity=Xjnorm(:,5:6);

%%types
Type0 = ones(size(Xi,1),1); % Type of student:
Type0((Xi(:,1) == 1) & (Xi(:,2) == 0)) = 2;
Type0((Xi(:,1) == 0) & (Xi(:,2) == 1)) = 3;
Type0((Xi(:,1) == 0) & (Xi(:,2) == 0)) = 4;

%% matrices for non-individual-specific terms
Inst = Data.TypeIndex;
MajorIndex=Data.TypeIndex(:,2);

X_MajorIndex = zeros(MCMC.J,length(unique(MajorIndex)));
for jj=1:MCMC.J
    X_MajorIndex(jj,MajorIndex(jj)) = 1;
end
assert(all( sum(X_MajorIndex') == 1));
X_MajorIndex = X_MajorIndex(:,2:end); %drop a major

Inst_FE_ind = 3;
X_InstitutionIndex = zeros(MCMC.J,length(unique(Inst(:,Inst_FE_ind))));
for jj=1:MCMC.J
    X_InstitutionIndex(jj,Inst(jj,Inst_FE_ind)) = 1;
end
assert(all( sum(X_InstitutionIndex') == 1));
inds_keep = find(~ismember(Data.InstName,'UNAB')); %drop Andres Bello dummy
X_InstitutionIndex = X_InstitutionIndex(:,inds_keep);


X_fixed = [XjArancel XjArancel.^2 XjSelectivity X_InstitutionIndex X_MajorIndex ~MCMC.platform];


%% match terms
X_interact_i = Sinorm(:,1:2);
X_interact_j = Xjnorm(:,1:2);
X_price_ij = XjArancel' - AidAmount(:,1);
X_price_ij(X_price_ij < 0) = 0; 
X_BVPdiscount_ij = zeros(size(X_price_ij));
X_BVPdiscount_ij(Over600 == 1, EverBVPj==1) = X_price_ij(Over600 == 1, EverBVPj==1);
X_dist_ij = Loci == Locj';


%% indiv-specific oo terms
typ_mat = (Type0 == [1 2 3]);
xloc = ((Loci == 13) | (Loci == 5) | (Loci == 8));  % Index 5
current_cohort = Data.current_cohort; %indicator for this year's hs class
is2011 = (Data.year==2011)*ones(size(xloc));
is2012 = (Data.year==2012)*ones(size(xloc));
X_oo = [ones(size(xloc)) Sinorm(:,1:2) xloc current_cohort is2011 is2012 AidAmount(:,1)]; 

%% terms to be multiplied by random coefficients
X_rc = Xjnorm(:,1:2);

%%
end